Cluster differential expressed genes

library_load <- suppressMessages(
    
    list(
        
        # Seurat 
        library(Seurat), 
        
        # Data 
        library(tidyverse), 
        library(openxlsx), 
        
        # Plotting 
        library(patchwork)
        
    )
)

# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")

# Source files
source("plotting_global.R")
source("bin/SeuratDEG.R")

# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()

options(warn=-1)

Parameter settings

so_file_nacl <- "data/object/seurat_sct_nacl.rds"
so_file_cpg <- "data/object/seurat_sct_cpg.rds"

Import Seurat

so_nacl <- readRDS(so_file_nacl)
so_cpg <- readRDS(so_file_cpg)

so <- list(so_nacl, so_cpg)

names(so) <- c("NaCl", "CpG")

Cluster annotation

ref <- readRDS("data/haemosphere/se_haemosphere.rds")

source("bin/SingleRQC.R")
so <- lapply(so, SingleRSeurat, ref=ref, cluster="SCT_snn_res.0.8")

Running SingleR in cluster mode

Running SingleR in cluster mode

set_factor_level <- function(so) {
    
    so$main_labels_cluster <- factor(so$main_labels_cluster, factor(names(color$main_labels_haemosphere)))
    so$fine_labels_cluster <- factor(so$fine_labels_cluster, factor(names(color$fine_labels_haemosphere)))
    
    return(so)
    
}

so <- lapply(so, set_factor_level)

Compute cluster differential expressed genes

# cluster_deg <- lapply(so, FindAllMarkers)

so_nacl <- so[["NaCl"]]
cluster_deg_nacl <- readRDS("result/cluster_deg/cluster_deg_nacl.rds")
# cluster_deg_nacl <- cluster_deg[["NaCl"]]
# saveRDS(cluster_deg_nacl, "result/cluster_deg/cluster_deg_nacl.rds")

so_cpg <- so[["CpG"]]
cluster_deg_cpg <- readRDS("result/cluster_deg/cluster_deg_cpg.rds")
# cluster_deg_cpg <- cluster_deg[["CpG"]]
# saveRDS(cluster_deg_cpg, "result/cluster_deg/cluster_deg_cpg.rds")

cluster_deg_nacl <- dplyr::left_join(cluster_deg_nacl, unique(so_nacl@meta.data[, c("SCT_snn_res.0.8", "main_labels_cluster", "fine_labels_cluster")]), by=c("cluster"="SCT_snn_res.0.8"))

cluster_deg_cpg <- dplyr::left_join(cluster_deg_cpg, unique(so_cpg@meta.data[, c("SCT_snn_res.0.8", "main_labels_cluster", "fine_labels_cluster")]), by=c("cluster"="SCT_snn_res.0.8"))

cluster_deg_nacl <- split(cluster_deg_nacl, f=cluster_deg_nacl$cluster)
cluster_deg_nacl <- cluster_deg_nacl[as.character(1:length(cluster_deg_nacl))]
names(cluster_deg_nacl) <- paste0("cluster_", names(cluster_deg_nacl))

cluster_deg_cpg <- split(cluster_deg_cpg, f=cluster_deg_cpg$cluster)
cluster_deg_cpg <- cluster_deg_cpg[as.character(1:length(cluster_deg_cpg))]
names(cluster_deg_cpg) <- paste0("cluster_", names(cluster_deg_cpg))

Export xlsx

# write.xlsx(cluster_deg_nacl, file="result/cluster_deg/cluster_deg_nacl.xlsx", overwrite=TRUE)

# write.xlsx(cluster_deg_cpg, file="result/cluster_deg/cluster_deg_cpg.xlsx", overwrite=TRUE)

Vulcaon plots

source("bin/SeuratQC.R")

vp_nacl <- lapply(cluster_deg_nacl, deg_volcano_plot, log2_thold=0.75, adjpvalue_thold=0.05)
vp_cpg <- lapply(cluster_deg_cpg, deg_volcano_plot, log2_thold=0.75, adjpvalue_thold=0.05)

NaCl

options(repr.plot.width=25, repr.plot.height=25)
wrap_plots(vp_nacl)

CpG

options(repr.plot.width=25, repr.plot.height=20)
wrap_plots(vp_cpg)

Heat maps

source("bin/SeuratDEG.R")

av_exp_nacl <- AverageExpression(so_nacl, assay="RNA", slot="data", group.by="SCT_snn_res.0.8")[["RNA"]]
av_exp_cpg <- AverageExpression(so_cpg, assay="RNA", slot="data", group.by="SCT_snn_res.0.8")[["RNA"]]

av_exp_nacl <- av_exp_nacl[, as.character(1:ncol(av_exp_nacl))]
av_exp_cpg <- av_exp_cpg[, as.character(1:ncol(av_exp_cpg))]

av_exp_hm_nacl <- lapply(cluster_deg_nacl, average_expression_hm, av_exp_nacl)
av_exp_hm_cpg <- lapply(cluster_deg_cpg, average_expression_hm, av_exp_cpg)

NaCl

options(repr.plot.width=25, repr.plot.height=25)
av_exp_grid_nacl <- gridExtra::arrangeGrob(grobs=av_exp_hm_nacl, ncol=5)
grid::grid.draw(av_exp_grid_nacl)

CpG

options(repr.plot.width=25, repr.plot.height=20)
av_exp_grid_cpg <- gridExtra::arrangeGrob(grobs=av_exp_hm_cpg, ncol=5)
grid::grid.draw(av_exp_grid_cpg)